Stochastic unraveling of time-local quantum master equations 
beyond the Lindblad class 
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A new method for stochastic unraveling of general time-local quantum master equations (QME) 
which involve the reduced density operator at time t only is proposed. The present kind of jump 
algorithm enables a numerically efficient treatment of QMEs which are not of Lindblad form. So 
it opens new large fields of application for stochastic methods. The unraveling can be achieved by 
allowing for trajectories with negative weight. We present results for the quantum Brownian motion 
and the Redfield QMEs as test examples. The algorithm can also unravel non-Markovian QMEs 
when they are in a time-local form like in the time-convolutionless formalism. 
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Quantum master equations (QMEs) are frequently 
used to describe time-independent as well as time- 
dependent phenomena in chemical physics, quantum op- 
tics, solid state physics, biological physics, etc. (see Ref. [I] 
for a number of typical examples). These QMEs de- 
scribe the time evolution of density matrices which are 
used in order to represent the mixed nature of the states. 
Stochastic unraveling is an efficient numerical tool for 
solving such equations. This method allows one to sim- 
ulate much larger and more complex systems with many 
degrees of freedom. It can, for example, be used to ac- 
curately describe femtochemical experiments in the liq- 
uid phase whose description has been limited until now 
to models with one or two effective interaction coordi- 
nates. In the unraveling scheme one considers an en- 
semble of stochastic Schrodinger equations (SSEs) which 
in the limit of a large ensemble resembles the respective 
QME. The numerical effort scales much more favorably 
with the size of the basis since one is now dealing with 
wave functions and not density matrices anymore (for a 
comparison of direct integrators, see Ref. ||). Another 
aspect of the stochastic methods is the possible physical 
interpretation of experiments detecting macroscopic fluc- 
tuations (e.g. photon counting) in various quantum sys- 
tems @. Most of the unraveling schemes §, |, |JJ| |, § 
have been restricted to QMEs of Lindblad form || which 
ensures that the reduced density matrix (RDM) stays 
positive semi-definite for all times and all parameters. 
Nevertheless there are many physical meaningful QMEs 
which result in positive- or almost positive-definite RDMs 
although they are not of Lindblad form. The increasing 
interest in descriptions beyond the Lindblad class such 
as the quantum Brownian motion [jlo| Ok the Redfield 
formalism non-Markovian schemes O, [T^ , fl5| , etc. 
resulted in various efforts to develop new stochastic wave 
function algorithms. 

Strunz et al. |l(], [if] developed the non-Markovian 
Quantum Diffusion Model. In general, this method 
can also be applied to QMEs in Markov approxima- 



tion even though they might not preserve positivity (see 
also Ref. |I(|). A similar approach was also proposed 
by Gaspard et al. Jf?]]. Very recently Stockburger and 
Grabert[^8) developed a method how to exactly repre- 
sent the RDM of a system coupled to a linear heat bath 
in terms of SSEs. The numerical properties of this ap- 
proach need to be explored. Breuer et al. [jl5| extended a 
scheme which they had used to calculate multi-time cor- 
relation functions to the unraveling of QMEs. Their 
technique is based on doubling the Hilbert space. In- 
stead of a single stochastic wave function one has a pair 
of them |]l5fl . This scheme conserves Hermiticity of the 
RDM only on average and not for every single realiza- 
tion. Thus, the deviation from Hermiticity is a quantity 
with statistical error and one has to perform a huge num- 
ber of realizations in order to achieve good convergence. 
Since stability and efficiency are crucial issues for unrav- 
eling algorithms we propose in this article an alternative 
approach which fulfills these criteria. 

The aim is to represent, in terms of quantum tra- 
jectories, the solution p(t) of a generalized time-local 
Hermiticity-conserving QME 



dp(t) 
dt 



= A{t)p{t) + P (t)A\t) 



M 



+ ^{C k (t)p(t)Et(t)+E k (t)p(t)Ct(t)}(l) 



k=l 



with the total number M of dissipative channels and ar- 
bitrary operators A(t), Ck{t), and Ek(t). Examples for 
these operatores are given below. Here we restrict the 
operators in such a way that the norm of the solution 
stays conserved. For readability we shall omit the time 
arguments in the following. 

In order to approach the problem let us define a state 
vector (IV'), \4>)) T spanning a doubled Hilbert space as 
proposed in Ref. [Hj. Unlike Ref. [t| the RDM shall be re- 
produced by an ensemble average (denoted by overbars) 



2 



of outer products of the vectors and \(f>): 



(2) 



A particular realization of the stochastic process will be 
denoted by the pair (\ip), <fi)). The averaging is per- 
formed over all trajectories possibly including a weighted 
sum over pure initial states. A vantage of this averaging 
is the conservation of Hcrmiticity for every single trajec- 
tory in contrast to Ref. |l5| . We note that this small mod- 
ification improves the numerical efficiency significantly. 

For the SSEs let us consider 2M independent possibly 
complex noise variables £L(i)- The superscripts denote 
which of the two terms from the Hermitian pair in the 
sum in Eq. (|l|) is taken and subscripts denote the various 
dissipative channels. All stochastic differentials ei££(i) 
are assumed to have zero mean, to be normalized and 
uncorrelated 120(1: 



<%i = 0, d£i*<%l = SijSkidt 



(3) 



Next, as an ansatz we construct a SSE which propagates 
the pair (|V>), \(f>)) 



M 2 



d\i>) = DM)dt + Y,Y, s ^) d & > ( 4a ) 

fe=l i=l 
M 2 

d\4>) = D 2 \4>)dt + Y,H s ^\4>) d ^ ■ ( 4b ) 



h=l i=l 



The operators Z?i and D 2 govern the deterministic and 
the operators Sj fe the stochastic part of the evolution. In 
general, they may depend on the state vector and explic- 
itly on time. After differentiating Eq. (Q), neglecting all 
terms higher than first order in dt, and assuming that 
ensemble averages always factorize Eli one obtains 



dp 



dt 



Di|VW|+I>2|MV| 

M 

Y, [slk\Wf\slt + slkWMst 



(5) 



fc=i 



dt + h.c. 



Comparing with Eq. (|l|) one notes that 5| fe has to equal 
Sf fc and S 2k has to equal Sl k . Moreover, one can see that 



S 2k - C k + 



a\ and 5| fc = E k + a k 



with a\ and a k being 

T 



arbitrary scalar functions of , \({>))~ L and possibly of 
time. Making the latter substitutions in Eq. (||) yields 
the constraint 



D 1= D 2 



A 



M 



fe=l 



a k *E k +ala 2 k *) . (6) 



Any quantum jump method is specified by jump rates 
pj, which have to be real scalar functions of (\ip), \(j>)). If 
n k {t) is the number of jumps in channel k and due to 
term i up to time t, the probability for n l k (t) to increase 



by one, i.e. the expectation value of both dn k and (dn k ) 2 , 
is equal to p k dt during the infinitesimal time interval dt. 
Thus, the noise variables obeying Eq. (||) are related 
"1 



to dn\(t) as 



"Sfc — 7=f e 



(7) 



The phase factor e %v does not change the RDM expres- 
sions within each realization and can be set to one. Sub- 
stituting Eq. (Q) into Eq. (|^) one finds that a k = — 
So the SSEs for our quantum jump method read 



d\iP) 




The jump rates p\ and p\ still remain free parameters. 
In the statistical limit their values have no influence on 
any averaged physical quantity. Nevertheless, it turns out 
that they can strongly influence the convergence behavior 
of the jump algorithm, i.e. they determine the statistical 
error of the observables calculated. A detailed discussion 
of this influence and utilization of such free parameters 
can be found in Ref. [2^. 

To ensure an efficient scheme with fast convergence one 
has to require that the norm of every single trajectory is 
constant in time. Asking for (4>\4>), ( l 0|V , )> e ^ c - being 
constant in time does not create a stable scheme but the 
condition of norm preservation of \i)}{4>\ + \4>) ("01 



(9) 



does. Unfortunately, applying this condition does not 
lead to positive values of the jump rates p l k for all tra- 
jectories at all times. However, since the p l k are arbi- 
trary real functions, they can be replaced by their ab- 
solute values. The price to pay is that we have to in- 
troduce an additional weight factor for the trajectories 
which jumps between one and minus one. In addition, 
there is a small deviation of the norm from unity be- 
cause in the regions where the p l k are replaced by their 
absolute values norm conservation is no longer guaran- 
teed. But in all our tests this deviation was far below 
one percent and neither effected numerical stability nor 
efficiency. The negative weights are actually needed to re- 
construct RDMs which are, in general, not positive-semi- 
definite. If the RDM stays positive-semi-definite during 
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its entire time-evolution the negative weights of some tra- 
jectories are not needed, i.e. all trajectories can be nor- 
malized to unity and represent physically pure states of 
the open quantum system. In the examples below the 
RDM can exhibit negative populations. This unphysical 
situation could probably be cured by applying an initial 
slippage to the initial state |23|, |24|. We note that this 
physically unreasonable RDMs occur because of unphys- 
ical initial states or because the QME is not physically 
correct or is applied in a parameter region where it is is 
not valid. Nevertheless an unraveling scheme has to be 
able to mimic also this unphysical behavior of the QME 
because in the ensemble average both should fully coin- 
cide. 

The condition applied to the QME (|l|) results in 
the additional constraint 



M 

£ 

k=l 



E\C k 



C\E k 



(10) 



and if applied to the deterministic part of the correspond- 
ing SSE (ph yields the total jump rate 



P 



\A + Ai\iP) + (il>\A + A^ 
<0|V^> + (V#) 



(11) 



All partial jump rates can be found subsequently making 
use of Eqs. (fTob and (O): 



Pi 



\cjE k \j>) + {mlc k \d>) 

\EjC k \^) + (MCjE^) 



(12a) 
(12b) 



In the rest of this article let us briefly show how the 
proposed method can be applied to two typical physical 
problems: the quantum Brownian motion and dissipative 
electron transfer within Redfield theory. In both cases 
the systems are described by Markovian QMEs which do 
not have Lindblad structure. The model of Brownian 
motion jjj describes a particle with mass to, coordinate 
g, momentum p and Hamiltonian Hg interacting with a 
thermal bath. In the high temperature limit of a bath of 
harmonic oscillators the relevant QME has the form 



dp 
~db 



[H S , P ] 



^7 

2h 



[<7> {P, P}] 



rwykT 



[q, hp}} (13) 



where 7 is the damping rate. Comparing with Eq. (|l]) 
one finds the operators of the jump algorithm (M = 2) 
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(14b) 
(14c) 
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FIG. 1: Time evolution of the third excited state of the 
harmonic oscillator in the quantum Brownian oscillator model 
for 7 = 10~ 3 lo, kT = 4.5w. The direct integration of the QME 
(thick solid line) is compared to the results of the quantum 
jump method with 1 trajectory (dot-dashed line), average of 
100 (thin solid line) and 1000 (broken line) trajectories. 



Modeling the particle as a harmonic oscillator with eigen- 
frequency u> one can compute the population dynamics 
depicted in Fig. |l|. The initial state of the oscillator is 
the pure state P33 = 1. As can be seen, the agreement of 
the results using our new stochastic method and using a 
direct integration of the QME is already quite good for 
one thousand samples. 

As a next test for the present quantum jump method 
we shall demonstrate the stochastic unraveling of the 
Redfield QME O, |§ 



[# s ,p] + -U [Ap,K } + [K,ptf] } 



(15) 



in which A is the relaxation operator and K the system 
part of the system-bath interaction |l^| . Let us consider 
a model for electron transfer in which the system includes 
a single reaction coordinate with the Hamiltonian |l2|, |2j| 

H s = #i|l)(l| +ff 2 |2)(2| +t>i 2 (|l)<2| + |2)(1|) (16) 



where H\ and H2 are the Hamiltonians of two cou- 
pled harmonic oscillators with frequency u>. We choose 
a potential configuration in the normal region with no 
barrier between the two harmonic potentials (change of 
free energy AE — 2oj, reorganization energy A = 3lu) 
with inter-center coupling V12 = u>. The bath is de- 
scribed by a cut-off frequency lo c = lo and temperature 
kT = uj/4. The system-bath interaction is characterized 
by the damping rate T = irr]/(A4 exp(l)) = cj/10 (see 
Rcf. [^U for details). After rearrangement of Eq. (|l5) 
one can easily identify the operators involved in Eq. (p ) 
(M = 1): C x = K, Ei = A, A = -iH s - KA. A 
Gaussian wave packet located at the donor state |1) and 
having energy slightly above the crossing of the har- 
monic potentials was chosen as initial state. The nu- 
merical simulation for about 1000 trajectories provides 
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COf/(27l) 



of every single trajectory. In this sense the jump rates 
were used as parameters to influence the efficiency. Neg- 
ative values for the weight of single trajectories allow 
for the reconstruction of non positive-semidefinite RDMs 
if required. The method was successfully tested for a 
simple electron transfer model and for Brownian motion 
and should allow for better quantum dynamical simula- 
tion of large systems. It can also unravel non-Markovian 
QMEs when they are in a time-local form like in the 
time-convolutionlcss formalism ]l3| or in methods using 
auxiliary density matrices to include the memory effects 
JUj as well as post-Markov master equations p3| . 



FIG. 2: Relaxation of the donor population for the electron 
transfer model. The solid line shows the exact solution of the 
QME, the dashed line one arbitrary trajectory, the dotted line 
an average over 500 trajectories. 
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FIG. 3: Occurrence of the expectation values of the pop- 
ulation on the donor state produced by the new unraveling 
scheme for the Redfield QME (dotted line) and the standard 
normalized jump method for the Lindblad QME (solid line) 
at time ut/{2n) — 3, both with 5000 trajectories. 



sufficiently converged and accurate results. Fig. || shows 
the relaxation of the ensemble averaged donor population 
Pi = (ip\l) (l\(f>) + ((f)\l) (l\ip) . A widely discussed prop- 
erty of the Redfield equation is that it does not conserve 
positivity Although Pi is always positive the tiny 

negative fraction in Fig. ^ is evidence for the existence 
of single realizations with negative Pi . The simulation of 
the same system within the so-called diabatic-damping 
approximation (25[ |2(| with a Lindblad QME by means 
of the standard quantum jump method | | ||, § 
keeps all values of Pi well confined between and 1. 

To summarize, a new method of stochastic unravel- 
ing of QMEs beyond the Lindblad form is proposed and 
thus new large fields of application for stochastic methods 
are opened. This progress became possible with the use 
of the wave-function pair in the doubled Hilbert space 
and the derivation of stable, almost normalized SSEs. 
The efficiency is determined by the behavior of the norm 
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